Destabilisation of the Subpolar North Atlantic prior to the Little Ice Age

The cooling transition into the Little Ice Age was the last notable shift in the climate system prior to anthropogenic global warming. It is hypothesised that sea-ice to ocean feedbacks sustained an initial cooling into the Little Ice Age by weakening the subpolar gyre circulation; a system that has been proposed to exhibit bistability. Empirical evidence for bistability within this transition has however been lacking. Using statistical indicators of resilience in three annually-resolved bivalve proxy records from the North Icelandic shelf, we show that the subpolar North Atlantic climate system destabilised during two episodes prior to the Little Ice Age. This loss of resilience indicates reduced attraction to one stable state, and a system vulnerable to an abrupt transition. The two episodes preceded wider subpolar North Atlantic change, consistent with subpolar gyre destabilisation and the approach of a tipping point, potentially heralding the transition to Little Ice Age conditions.


Supplementary Figure 2. Trends in AR(1) after detrending using a short bandwidth.
Records are detrended using a bandwidth of 12 years to remove the irregular oscillations that characterise each episode. The annual values and the estimated trend (red line) are represented on the upper row. The middle row shows the residuals after detrending and the bottom row the AR(1) values obtained using a window length of 70 years. a) δ 18 O and b) Growth Index records during the first episode. c) δ 18 Oshell during the second episode. Declining resilience is detected over the yellow and blue shaded intervals corresponding to the first and second episodes respectively. The hatched region comprises the initial window interval for each episode.

Supplementary Figure 3. a)
Reconstructed sea bottom temperature from benthic foraminiferal δ 18 O at Malangen Fjord as a proxy for the influence of the Atlantic warm inflow 1 . b) Sortable silt mean grain size (SS ) record from South Iceland as a proxy for near-bottom flow speed of the Iceland-Scotland Overflow Water (ISOW), which is influenced by the strength of the warm Atlantic inflow 2 . c) Difference in δ 18 O measured on the foraminiferal species N. pachyderma and T. quinqueloba (Δδ 18 ONps-Tq) in the Labrador Sea as a proxy for the influence of warm Atlantic waters influenced by the SPG strength 3 . d) Location of the proxy records shown on the left figure (a-c) and regional ocean circulation. The colour map represents April sea surface temperatures during 2016, obtained from the ESA Climate Change Initiative data 4 . The land image was obtained from the NOAA National Oceanographic Data Center 5 .

Effect of the variable number of replicates and ontogenetic trends on the resilience trends
Before estimating changes in resilience, we reviewed the methods used to construct each proxy to assess their potential effects in autocorrelation and variance. The isotope series were obtained by averaging the number of replicate samples each year, which varies between one and six throughout both records 15,16 . In this case, we could expect intervals with a smaller sample size to have the largest variance. In both series, the number of replicates before 1170 CE is mostly one and increases to two between 1170 and 1240 CE 15,16 . If the sample size effect were the primary factor driving the variance trend, it would decrease over the first episode as the number of replicates increases; instead, a positive trend is observed.
The construction of shell-growth records involves more steps because shell-growth is influenced by age: young bivalves grow faster and produce wider and more variable annual increment bands compared with those produced during the mature years. Butler et al. 17 built two versions of the same record using different approaches to remove the age-related trend. The first method fits a negative exponential curve to the measurements of each shell after stabilising variance, whereas the second one follows the Regional Curve Standardisation (RCS) method developed in dendrochronology. Both versions exhibit similar trends in AR(1) and variance ( Supplementary Fig. 4c,d), indicating that neither approach affects the trends in the indicators. However, the variance stabilisation methods do not entirely remove the ontogenetic trend. We assess the persistence of the age-related trend by computing the resilience metrics over the individual shell-growth records. Before computing the trends in the indicators, the agerelated trend is removed by fitting a negative exponential curve after stabilising variance using a data-adaptive power-transformation method 18 . The results indicate that the ontogenetic trend in variance persists, resulting in significant negative trends in this indicator ( Supplementary Fig. 5b). The earliest years are also slightly more autocorrelated, resulting in small negative trends in AR(1) (Supplementary Fig. 5a). In this scenario, the age-related trends might modify the trends in the indicators over the intervals where new shell records are introduced to the chronology. The first episode coincides with the successive introduction of five shells ( Supplementary Fig. 4e), possibly altering the trend in the indicators. To investigate whether the trends respond to the introduction of new shells, we analysed the individual shell series around 1200 CE. The individual series share a common trend and similar values for both metrics between 1200 and 1260, regardless of their age ( Supplementary Fig. 6), suggesting that the ontogenetic trends do not drive the observed signal.
In addition to the age-related effects, the variance trends can also be affected by the variable number of replicates through time, expecting to observe increased variance over regions with the smallest sample size. To assess the probability of obtaining the observed signal due to the combination of both effects, we compared the observed trends to those obtained from three null models. Each null model comprises 1,000 surrogate series, and each series was created by introducing randomly chosen individual shell records at the same time as in the original chronology between 1100 and 1260. The probability of obtaining the observed trend under each scenario is measured as the proportion of surrogate series that yield higher or equal Kendall τ values as those obtained from a series built in the same way but preserving the original order in the shell records.
The first null model aims to test only the age-related effects without considering the effect of the variable number of replicates. In this model, the indicators are computed over the individual random records and annually-averaged to obtain each surrogate series. The results from this experiment indicate that it is possible to obtain positive trends in AR(1) and variance purely associated with age-related effects ( Supplementary  Fig. 7a). However, the proportion of series that yield larger values than the measured in the original records is less than 0.02 for both metrics.
The second experiment assesses the effect of the variable number of replicates through time on the indicators. In this case, each surrogate series is generated by averaging synthetic white-noise series of the same length, mean and variance as the shell records from the original chronology. This experiment demonstrates that the variable number of replicates significantly affects the variance trends ( Supplementary Fig. 7b), resulting in negative trends because the number of replicates starts with one and increases progressively to six before 1260 ( Supplementary Fig. 4e).
The third null model tests the effect of both factors combined. Each surrogate series is created by averaging the individual records using a biweight robust mean, as in the original records. The resilience indicators are then computed over the resulting series. In this experiment, the trends in AR(1) do not exhibit a positive bias, suggesting that autocorrelation is destroyed by averaging the series unless the individual records contain the same signal. On the other hand, variance is slightly skewed towards negative trends ( Supplementary Fig. 7c), indicating that the effect of the variable number of replicates outweighs the age-related effects. The proportion of surrogate series that yield higher values than the observed in the original record is less than 0.002 for each metric. The results from the null models indicate that it is improbable that the observed trend during the first episode resulted from age-related effects or due to the variable number of replicates, supporting the idea that the environment is the primary driver of the observed signal before 1260. c) The individual shell-growth records are averaged to obtain each surrogate series before computing the resilience metrics. The red dotted vertical line indicates the Kendall τ value measured on a series built following the process for each experiment but preserving the original order of the individual shell-growth records.